
require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("psych")
require(psych)
library(gridExtra)
library(grid)

# Please make sure to change the path to where the appropriate dataset is located on your computer

sc.val <-  read.csv("PNAS_VAL.csv") 

# Columns

sponsor.require <- 1
p.tourist.visa <- 2
p.fakedocs <-  3
k <- 4
my.seed <- 5

legaltourist.attempts <- c(6:25)
ability.legal.cols <- c(26:86)
ability.fakedocs.cols <- c(87:147)
ability.work.cols <- c(148:208)

tot.mig <- c(209:269)
legal.mig <- c(270:330)
fd.mig <- c(331:391) # illegal
ut.mig <- c(392:452) # semi-legal
stud.mig <- c(453:513)
hs.mig <- c(514:574)
ls.mig <- c(575:635)
fam.mig <- c(636:696)

legal.start <- 232
ut.start <- 9
fd.start <-  9
tot.start <-  250

####

table(sc.val[,k])

### Function for shared legend

grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = grid::unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           heights = grid::unit.c(unit(1, "npc") - lheight, lheight)))
  grid.newpage()
  grid.draw(combined)
  
}


### Visa area plot

cat <- c(rep("Student",12), rep("High-Skilled",12),rep("Low-Skilled",12), rep("Family",12))

stud.means <- t(as.vector(colMeans(sc.val[which(sc.val[,k] == 1),stud.mig[1:12]])))
hs.means <- t(as.vector(colMeans(sc.val[which(sc.val[,k] == 1),hs.mig[1:12]])))
ls.means <- t(as.vector(colMeans(sc.val[which(sc.val[,k] == 1),ls.mig[1:12]])))
fam.means <- t(as.vector(colMeans(sc.val[which(sc.val[,k] == 1),fam.mig[1:12]])))

cat <-  c(rep("Student", 12), rep("High-Skilled",12), rep("Low-Skilled",12), rep("Family",12))
years <- c(rep(1:12,4))
Value1 <- append(stud.means,hs.means)
Value2 <- append(ls.means,fam.means)
all <- append(Value1,Value2)
v = seq(from = 1, to =12, by = 1)

data <-  data.frame(cat,years,all)

data <- data[order(data$years),]

r <- ggplot(data, aes(x=data$years, y=data$all, group=cat, fill=cat)) + 
  geom_area(position="fill") + xlab("Years") + ylab("Proportion of Migrants Granted Visa") + theme_minimal() +theme(legend.title = element_blank()) +
  scale_x_continuous(name="Years", breaks = v) + scale_fill_brewer(palette="Set3") + ggtitle("(a) ABM")

# Data from DHS

stud.dhs <- c(8274, 8073, 7407, 6856, 7311, 8831, 6236, 11329, 9798, 9037, 8388, 7267)
hs.dhs <-  c(1006, 1275, 1439, 1780, 1823, 2179, 2151, 2063, 1806, 1748, 1533, 1247)
ls.dhs <-  c(12821, 12150, 13042, 11315, 11943, 14864, 18918, 12896, 7849, 7624, 8745, 8376)
fam.dhs <-  c(6141, 4803, 5696, 5211, 5032, 6218, 4988, 5790, 4833, 5386, 6074, 5818)

dhs.1 <-  append(stud.dhs,hs.dhs)
dhs.2 <- append(ls.dhs,fam.dhs)
dhs.all <-  append(dhs.1, dhs.2)

dhs.data <-  data.frame(cat, years,dhs.all)
dhs.data <- dhs.data[order(dhs.data$years),]

v = seq(from = 1, to =12, by = 1)

s <- ggplot(dhs.data, aes(x=dhs.data$years, y=dhs.data$dhs.all, group=dhs.data$cat, fill=dhs.data$cat)) + 
  geom_area(position="fill")+ xlab("Years") + ylab("") + theme_minimal() + theme(legend.title = element_blank()) +
  scale_x_continuous(name="Years", breaks = v) + scale_fill_brewer(palette="Set3") + ggtitle("(b) DHS")

pdf("FigC2", width = 10, height = 7, onefile = FALSE)
grid_arrange_shared_legend(r,s,nrow=1,ncol=2)
dev.off()

